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ANGIOGRAPHY METHOD AND APPARATUS 

Background of the Invention 

The present invention relates to the medical 
imaging arts. It particularly relates to angiography 
using the magnetic resonance imaging (MRI) and computed 
tomography (CT) medical imaging techniques, and will be 
described with particular reference thereto. However, the 
invention will also find application in conjunction with 
other three-dimensional imaging modalities as well as in 
other imaging arts in which thin structures or networks 
with overlapping or furcated portions are advantageously 
differentiated from extraneous imaged structures and 
background noise or tracked in three dimensions. 

Catastrophic medical events such as heart 
attacks and strokes that result from underlying vascular 
problems are a leading cause of death in the United 
States. Plaque buildup on the inside of the vascular 
walls can lead to strokes, coronary heart disease, and 
other medical conditions such as vessel rupture. Many 
Americans suffer from chronic vascular diseases which 
degrade quality of life. 

Angiography relates to the imaging of blood 
vessels and blood vessel systems, and as such is a 
powerful medical diagnostic for identifying and tracking 
vascular diseases. Angiography enables improved surgical 
planning and treatment, improved diagnosis and convenient 
non-invasive monitoring of chronic vascular diseases, and 
can provide an early indication of potentially fatal 
conditions such as aneurysms and blood clots. The ability 
of certain types of angiography to accurately characterize 
the vessel lumen is particularly valuable for diagnosing 
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plaque buildup on the vessel walls. 

Angiography is performed using a number of 
different medical imaging modalities, including biplane 
X-ray/DSA, magnetic resonance (MR) , computed tomography 
5 (CT) , ultrasound, and various combinations of these 
techniques. Two-dimensional or three-dimensional 

angiographic data can be acquired depending upon the 
medical imaging modality and the selected operating 
parameters. Certain types of angiography employ invasive 

10 contrast enhanced methodologies in which a contrast agent 
that accentuates the vascular image contrast is 
administered to the patient prior to the imaging session. 
Some angiography techniques, such as MR imaging, are also 
capable of providing vascular contrast using non-invasive 

15 methodologies that take advantage of intrinsic aspects of 
the vascular system, such as the blood motion or flow, to 
enhance the vascular contrast without employing an 
administered contrast agent. For either contrast-enhanced 
or non-contrast-enhanced angiography, the vasculature 

2 0 imaging is effectuated by either a signal enhancement in 
the vascular regions (white blood angiography) , or by a 
signal suppression in the vascular regions (black blood 
angiography) . 

The analysis of angiographic images by medical 

2 5 personnel is often hindered by image imperfections or 

intervening non-vascular structures (e.g., bone, organs, 
and the like) . Even in the absence of such problems, 
however, the sheer complexity of the vascular system and 
its myriad sub-systems severely complicates image 

3 0 interpretation. 

With reference to FIGURE 1, a schematic portion 
of an exemplary vasculature is shown, including an 
arterial sub-system A and a venous sub-system V. As is 
often the actual situation in the human body, the two sub- 
35 systems A, V are shown in FIGURE 1 arranged in a 
substantially parallel manner. Furthermore, there are 
numerous points where, in the view shown, an artery 
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portion overlaps a vein portion: exemplary points are 
designated AV. Similarly, there are numerous points where 
the a vein portion overlaps an artery portion: exemplary 
points are designated VA. Another complexity arises at 
5 furcation points. FIGURE 1 shows exemplary artery 
bifurcations AB and exemplary vein bifurcations VB. 

With reference to FIGURE 2, an exemplary 
vascular crossing is shown, in which a vessel V 1 and a 
vessel V 2 cross. In three-dimensional angiography, the 

10 image is typically created by imaging a plurality of 
parallel planes S which are then combined to form a three- 
dimensional image representation. For an exemplary slice 
S o oriented perpendicular to the vessels V, and V 2 , the 
image of the vessel V, in the plane S Q is shown 

15 superimposed as W v Similarly the image of the vessel V 2 
in the plane S Q is shown superimposed as W 2 . 

With reference to FIGURE 3A, it is seen that the 
vessel images W, and W 2 are overlapping in the image slice 
S o . FIGURE 3B shows the overlapping vessel images W 1 and 

20 W 2 as they would appear in a typical angiographic image. 
Since the contrast is essentially identical for W 1 and W 2 , 
it will not be clear to medical personnel whether the 
overlapping vessels represent crossing arteries, crossing 
veins, an artery crossing a vein, a vein crossing an 

25 artery, a vein bifurcation point, an artery bifurcation 
point, a closely positioned but non-overlapping pair of 
vessels, et cetera. 

The vascular complexity issues described with 
reference to FIGURES 1 through 3B are advantageously 

3 0 addressed by automated vascular segmentation systems. 
These systems differentiate the vasculature from non- 
vascular structures, background levels, imaging system 
artifacts such as noise, and the like. Many segmentation 
engines employ tracking systems which track a vessel 

3 5 starting from an initial seed location. Tracking systems 
can track the vessel skeleton while simultaneously 
quantifying the vessel lumen, and such systems are 


particularly useful for accommodating the varying vessel 
diameters usually encountered in following a blood vessel. 
Tracking systems also can separate out individual vessel 
branches. In the exemplary FIGURE 1, a tracking system 
5 starting at artery seed AS will track the arterial branch 
A, while the tracking system starting at vein seed VS will 
track the venous branch V. In this manner, artery-vein 
separation is achievable. 

However, tracking methods of the prior art have 

10 numerous disadvantages, principally due to the localized 
nature of the tracking analysis. Tracking can be 
cumbersome and inaccurate, particularly in areas of very 
high vascular densities such as in the brain. Bifurcation 
points, tortuous or occluded vessels, vessel overlaps, 

15 intertwined vessels, partial volume averaging and other 
imaging artifacts, and vessel gaps can act alone or in 
various combinations to produce localized disjoints of 
the vascular path (or localized disjoints of the 
angiographic image of the vascular path) which prevent 

20 successful vessel tracking. Furthermore, at vessel 
overlaps the wrong vascular system may be tracked. For 
example, a tracking system following the arterial branch 
A of FIGURE 1 could fail and begin tracking the venous 
branch V at any of the crossing points AV, VA. 

2 5 The present invention contemplates an improved 

angiographic method and apparatus which overcomes the 
aforementioned limitations and others. 

summary of the Invention 

According to one aspect of the invention, an 

3 0 apparatus is disclosed for producing an angiographic image 

representation of a subject. An imaging scanner acquires 
imaging data from at least a portion of a subject. The 
imaging data includes vascular contrast. A reconstruction 
processor reconstructs an image representation from the 
35 imaging data. The image representation is formed of image 
elements and exhibits vascular contrast. A processor 
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converts the image representation into an edge-enhanced 
image representation having enhanced vascular edges and 
divides the edge-enhanced image representation into at 
least one two-dimensional slice formed of pixels. For 
5 each slice, the processor flood-fills the vascular edges 
to form filled regions defined by pixels having a first 
value, identifies vessel centers through iterative removal 
of pixels having the first value from around the edges of 
the filled regions, and segments, tracks, extracts, 

10 enhances, or identifies vascular information contained in 
the angiographic image using the identified vessel centers 
as operative inputs. 

According to another aspect of the invention, a 
method is disclosed for characterizing a vascular system 

15 in a three-dimensional angiographic image comprised of 
voxels. A two-dimensional slice formed of pixels is 
extracted from the angiographic image. Imaged vascular 
structures are located in the slice. The imaged vascular 
structures are flood-filled to form filled regions defined 

2 0 by pixels having a first value. The edges of the filled 
regions are iteratively eroded to identify vessel center 
points. The extracting, locating, flood-filling, and 
eroding are repeated for a plurality of slices to generate 
a plurality of vessel center points that are 

2 5 representative of the vascular system. 

According to another aspect of the invention, a 
method is disclosed for tracking a vascular system in an 
angiographic image. A plurality of vessel centers are 
identified in three dimensions that are representative of 

30 the vascular system. A first vessel center is selected. 
A first vessel direction is found corresponding to the 
local direction of the vessel at the first vessel center. 
A first slice is defined that is orthogonal to the first 
vessel direction and includes the first vessel center. 

35 Vessel boundaries are estimated in the first slice by 
iteratively propagating a closed geometric contour 
arranged about the first vessel center. The selecting, 
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finding, defining, and estimating are repeated for the 
plurality of vessel centers. The estimated vessel 
boundaries are interpolated to form a vascular tree. 

According to yet another aspect of the 
5 invention, an apparatus for characterizing a vascular 
system in a three-dimensional angiographic image comprised 
of voxels is disclosed. A means is provided for 
extracting from the angiographic image a two-dimensional 
slice formed of pixels. A means is provided for locating 

10 imaged vascular structures in the slice. A means is 
provided for flood-filling the imaged vascular structures 
to form filled regions defined by pixels having a first 
value. A means is provided for iteratively eroding the 
edges of the filled regions to identify vessel center 

15 points. A means is provided for generating a plurality of 
vessel center points that are representative of the 
vascular system. The means for generating is in operative 
communication with the means for extracting, the means for 
locating, the means for flood-filling, and the means for 

2 0 eroding. 

According to still yet another aspect of the 
invention, an apparatus for tracking a vascular system in 
an angiographic image is disclosed. A means is provided 
for identifying a plurality of vessel centers in three 

2 5 dimensions that are representative of the vascular system. 

A means is provided for selecting a first vessel center. 
A means is provided for finding a first vessel direction 
corresponding to the local direction of the vessel at the 
first vessel center. A means is provided for defining a 
30 first slice that is orthogonal to the first vessel 
direction and includes the first vessel center. A means 
is provided for estimating vessel boundaries in the first 
slice by iteratively propagating a closed geometric 
contour arranged about the first vessel center. A means 

3 5 is provided for interpolating the estimated vessel 

boundaries to form a vascular tree after the selecting, 
finding, defining, and estimating have been repeated for 


- 7 - 

the plurality of vessel centers. 

One advantage of the present invention is that 
it is a global technique which overcomes segmentation 
failures often encountered at vessel disjoints and 
5 overlaps by localized techniques. 

Another advantage of the present invention is 
that it is compatible with both white blood angiography 
and black blood angiography. 

Another advantage of the present invention is 
10 that it provides rapid and accurate vessel boundary 
estimation using propagation of geometric contours. The 
propagating can advantageously incorporate the gray scale 
image information through fuzzy membership classification 
of image elements in the neighborhood of the contour. 
15 Another advantage of the present invention is 

that it provides locations of furcation and vessel overlap 
points globally throughout the angiographic volume. This 
information can be used by vessel trackers or other vessel 
segmentation systems to improve accuracy and speed. 
2 0 Yet another advantage of the present invention 

is that it improves tracking speed and accuracy by 
providing a plurality of vessel center points or tags that 
are representative of the global vascular system. 

Still yet another advantage of the present 

2 5 invention is that it advantageously retains the sharp 

vascular edges and accurate vessel lumen information 
typically achieved by black blood data acquisition during 
the vessel segmentation processing. 

Numerous additional advantages and benefits of 

3 0 the present invention will become apparent to those of 

ordinary skill in the art upon reading the following 
detailed description of the preferred embodiment. 

Brief Description of the Drawings 

The invention may take form in various 
35 components and arrangements of components, and in various 
steps and arrangements of steps. The drawings are only 
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for the purpose of illustrating preferred embodiments and 
are not to be construed as limiting the invention. 

FIGURE 1 shows a schematic drawing of a pair of 
intertwined vascular structures that have several 
5 bifurcation points; 

FIGURE 2 shows a schematic drawing of a vessel 
crossing with typical angiographic slice images 
superimposed ; 

FIGURE 3A shows a typical angiographic slice in 
10 the overlap region of FIGURE 2 with the vessels 
artificially differentiated; 

FIGURE 3B shows a typical angiographic slice in 
the overlap region of FIGURE 2 without the artificial 
differentiation ; 
15 FIGURE 4 shows an exemplary magnetic resonance 

angiography (MRA) system formed in accordance with an 
embodiment of the invention; 

FIGURE 5 shows an overview of an exemplary 
angiographic segmentation method formed in accordance with 
2 0 an embodiment of the invention; 

FIGURE 6A schematically shows a plurality of 
vessel centers that are representative of a bifurcated 
vascular branch; 

FIGURE 6B schematically shows a plurality of 

2 5 estimated vessel boundaries corresponding to the vessel 

centers of FIGURE 6A; 

FIGURE 6C schematically shows a representation 
of the bifurcated vascular branch reconstructed using the 
vessel centers and boundaries of FIGURES 6A and 6B; 

3 0 FIGURE 7 shows an exemplary embodiment of the 

pre-processor of the segmenter of FIGURE 5; 

FIGURE 8 shows an exemplary embodiment of the 
pseudo-white blood angiographic (pseudo-WBA) processor of 
FIGURE 7; 

3 5 FIGURE 9 shows an exemplary embodiment of the 

bone-air-muscle mask processor of FIGURE 8; 

FIGURE 10 shows an exemplary embodiment of the 
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class assignment processor of FIGURE 9; 

FIGURE 11 shows an exemplary embodiment of the 
three-dimensional high curvature feature processor of 
FIGURE 5; 

5 FIGURE 12 shows an exemplary embodiment of the 

decomposition of the Gaussian derivative kernel into 
cosinusoidal and sinusoidal components; 

FIGURE 13 shows an exemplary time-of-f light 
white blood angiographic image of the carotid area that 
10 includes five isolated vessels as well as one pair of 
vessels emerging from a bifurcation; 

FIGURE 14 shows the edge volume corresponding to 
the image of FIGURE 13; 

FIGURES 15 and 16 show an exemplary embodiment 
15 of the vessel centers processor of FIGURE 5; 

FIGURE 17A schematically shows an exemplary 
vessel structure corresponding to a single vessel center; 

FIGURES 17B, 17C, and 17D schematically show the 
recursive eroding of the vessel structure of FIGURE 17A; 
2 0 FIGURE 18A schematically shows an exemplary 

vessel structure corresponding to a pair of overlapping 
vessel centers; 

FIGURES 18B, 18C, and 18D schematically show the 
recursive eroding of the vessel structure of FIGURE 18A; 

2 5 FIGURE 19 schematically shows the decomposition 

of a sguare moving window for use in recursive erosion 
into an upper left component and a lower right component; 

FIGURE 2 0 shows an exemplary method for the 
first pass of an exemplary two-pass recursive erosion 

3 0 process; 

FIGURE 21 shows an exemplary method for the 
second pass of an exemplary two-pass recursive erosion 
process ; 

FIGURE 22A shows a synthetic structure 
35 consisting of two isolated circular flood-filled regions; 

FIGURE 22B shows the results of the recursive 
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erosion method of FIGURES 2 0 and 21 as applied to the 
synthetic structure of FIGURE 2 2A; 

FIGURE 22C shows a synthetic structure 
consisting of two overlapping circular flood-filled 
5 regions; 

FIGURE 22D shows the results of the recursive 
erosion method of FIGURES 20 and 21 as applied to the 
overlapping synthetic structure of FIGURE 22C; 

FIGURE 2 3 shows an exemplary path tracking 
10 process formed in accordance with an embodiment of the 
invention that employs a representative set of vessel 
center tags; 

FIGURE 24 shows a robust general tracking system 
formed in accordance with an embodiment of the invention 
15 that employs a representative set of vessel center tags to 
ensure complete tracking of the full vascular skeleton; 

FIGURE 25 shows an exemplary embodiment of the 
direction processor of FIGURE 23; 

FIGURE 26A schematically shows the exemplary 
2 0 vessel structure of FIGURE 18A with initial geometric 
contours arranged about the vessel centers; 

FIGURE 2 6B schematically shows the propagating 
of the initial geometric contours of FIGURE 2 6A; 

FIGURE 27 shows an exemplary embodiment of the 
25 vessel boundaries contour fitting; 

FIGURE 28 shows the exemplary embodiment of the 
contour propagation processor of FIGURE 27; and 

FIGURE 29 shows an exemplary embodiment of a 
three-dimensional display processor for displaying the 
30 vascular skeleton generated by the method of FIGURE 24. 

Detailed Description of the Preferred Embodiments 

With reference to FIGURE 4 , a magnetic resonance 
imaging system that suitably practices angiographic 
imaging in accordance with an embodiment of the invention 
35 is described. Although the invention is described herein 
with respect to a magnetic resonance imaging embodiment, 
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those skilled in the art will appreciate that the 
invention is applicable to a broad range of angiographic 
modalities and techniques, including but not limited to 
contrast-enhanced magnetic resonance angiography, non- 
5 contrast enhanced magnetic resonance angiography, computed 
tomographic angiography, and fused magnetic 
resonance/computed tomography angiographic techniques. 
The invention is also suitably practiced in conjunction 
with either white blood angiography (WBA) or black blood 

10 angiography (BBA) . 

With reference to FIGURE 4, a magnetic resonance 
imaging (MRI) scanner 10 typically includes 
superconducting or resistive magnets 12 that create a 
substantially uniform, temporally constant main magnetic 

15 field B Q along a z-axis through an examination region 14. 
Although a bore-type magnet is illustrated in FIGURE 4, 
the present invention is equally applicable to open magnet 
systems and other known types of MRI scanners. The magnets 
12 are operated by a main magnetic field control 16. 

2 0 Imaging is conducted by executing a magnetic resonance 

(MR) sequence with the subject being imaged, e.g. a 
patient 42 in a magnetic resonance angiography (MRA) 
session, placed at least partially within the examination 
region 14, typically with the region of interest at the 
25 isocenter. 

The magnetic resonance sequence entails a series 
of RF and magnetic field gradient pulses that are applied 
to the subject to invert or excite magnetic spins, induce 
magnetic resonance, refocus magnetic resonance, manipulate 

3 0 magnetic resonance, spatially and otherwise encode the 

magnetic resonance, to saturate spins, and the like. More 
specifically, gradient pulse amplifiers 20 apply current 
pulses to a whole body gradient coil assembly 22 to create 
magnetic field gradients along x-, y- , and z-axes of the 
35 examination region 14. 

An RF transmitter 24, preferably digital, 
applies RF pulses or pulse packets to a whole-body RF coil 
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2 6 to transmit RF pulses into the examination region. A 
typical RF pulse is composed of a packet of immediately 
contiguous pulse segments of short duration which taken 
together with each other and any applied gradients achieve 
5 a selected magnetic resonance manipulation. The RF pulses 
are used to saturate, excite resonance, invert 
magnetization, refocus resonance, or manipulate resonance 
in selected portions of the examination region. 

For whole-body applications, the resulting 

10 resonance signals, generated as a result of a selected 
manipulation, are also picked up by the whole-body RF coil 
26. Alternately, for generating RF pulses in limited 
regions of the subject, local RF coils are placed 
contiguous to the selected region. For example, as is 

15 known in the art, an insertable head coil 2 8 is inserted 
surrounding a selected brain region at the isocenter of 
the bore. Other surface coils or other such specialized 
RF coils may also be employed. For example, the RF system 
optionally includes a phased array receive coil (not 

20 shown) which enables partial parallel imaging (PPI) 
techniques known to the art. In one suitable embodiment, 
the whole-body RF coil 2 6 induces resonance and the local 
RF coil or coil array receives magnetic resonance signals 
emanating from the selected region. In other embodiments, 

25 the local RF coil both excites and receives the resulting 
magnetic resonance signals. 

Regardless of the RF coil configuration and the 
application thereof, the resultant RF magnetic resonance 
signals that are picked up by one or another of the RF 

30 coils is received and demodulated by an RF receiver 32. 
A sequence control processor 34 controls the gradient 
pulse amplifiers 20, the RF transmitter 24, and the RF 
receiver 32 to produce integrated MRI pulse sequence and 
readout waveforms that generate the magnetic resonance 

3 5 (MR) signals and optional echoes, provide appropriate 
encoding gradients to spatially encode the resultant MR 
response, and coordinate MR pickup and receive operations. 
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The MR I sequence typically includes a complex 
series of magnetic field gradient pulses and/or sweeps 
generated by the gradient amplifiers 20 which along with 
selected RF pulses generated by RF coils 26, 28 result in 
5 magnetic resonance echoes that map into k-space. The 
resultant magnetic resonance data is stored in a k-space 
memory 36. The k-space data is processed by a 

reconstruction processor 38, which is typically an inverse 
Fourier transform processor or other reconstruction 

10 processor known to the art, to produce a reconstructed 
image representation that is stored in an image memory 40. 

In magnetic resonance angiography (MRA) , a 
patient 42 is imaged by the MRI system 10 using imaging 
conditions that provide either an enhanced high intensity 

15 signal from the vascular regions (white blood angiography) 
or a suppressed low intensity signal from the vascular 
regions (black blood angiography) for the vascular 
regions. The enhanced or suppressed vascular signal 
provides contrast for the vascular system in the resultant 

20 image. In the exemplary embodiment of FIGURE 4, the 
carotid area of the patient 42 is imaged. Optionally, the 
patient receives a magnetic resonance contrast agent 44 to 
improve the vascular contrast, i.e. contrast-enhanced MRA 
is performed. An angiographic sequence such as a time of 

25 flight (TOF) white blood sequence, a gradient recalled 
black blood sequence which can be used in combination with 
large bipolar gradients or pre-saturation RF pulses, a 
spin-echo (SE) type black blood sequence which utilizes 
wash-out effects and which does not require large bi-polar 

3 0 gradients or pre-saturation RF pulses, or other 
appropriate angiographic sequence is employed. The 
selected angiographic sequence is applied by the sequence 
control processor 34 and effectuates acquisition of the 
MRA data. 

3 5 Regardless of the choice of MRA imaging 

modality, the k-space data is collected in the k-space 
memory 3 6 and reconstructed by the reconstruction 
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processor 38 to generate the MRA volume image data 40. 
The MRA volume image data is in the form of a three- 
dimensional gray scale image representation of the 
examined area of the patient, which has good contrast for 
5 the vascular system relative to other body tissues of the 
patient 42. Typically a three dimensional image data 
comprising multiple slices is acquired to provide volume 
imaging of the patient 42. However, it is also 
contemplated that only a single image slice is acquired in 

10 a patient imaging session. 

A pre-processor 4 6 conforms the angiographic 
data to a standardized format preparatory to further 
processing. For example, the pre-processor 46 optionally 
smooths the data, converts the data into an isotropic 

15 format, re-sizes the image, performs pre-f iltering to 
remove noise or extraneous features such as non-vascular 
contrast that interferes with the vascular image, and the 
like. For BBA, the pre-processor preferably intensity- 
inverts the image so that the vascular areas are 

2 0 represented by high intensities. 

An edge volume processor 48 receives the output 
of the pre-processor 46, and applies a mathematical 
transformation such as second order spatial 
differentiation to obtain an edge-enhanced volume that 

2 5 particularly emphasizes the edges of the vasculature. 

A vessels center tagger 50 receives the output 
of the edge volume processor 48 and searches the edge 
volume for vessel centers, e.g. on a slice-by-slice 
basis. The located vessel center points in each plane are 

3 0 tagged. Vessel centers corresponding to overlapping 

vessel images are advantageously identified and 
particularly noted. The vessel centers along with the 
overlap tags are supplied along with the edge volume image 
and the pre-processed volume image to a vessel 
3 5 segmentation engine 52. 

The segmentation processor 52 processes the 
angiographic data to segment the vascular structure. The 
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segmenting can be performed in several ways, including: 
identifying the vasculature as a set of voxels 
corresponding to the imaged vascular structures; 
separating of the vascular regions of the image from the 
5 non-vascular regions of the image; assigning a first value 
to voxels corresponding to imaged vascular structures, and 
a second value to voxels corresponding to imaged non- 
vascular structures; and removing contrast corresponding 
to non-vascular structures from the image. 

10 The segmented vascular information is preferably 

graphically displayed on an appropriate user interface 54, 
typically as a two-dimensional or three-dimensional 
graphical image representation. The vascular information 
can of course also be stored electronically or on a 

15 magnetic storage medium, printed on paper, and et cetera 
(not shown) . 

With reference to FIGURE 5, a method 70 for 
post-acguisition processing of acguired raw angiographic 
volume image data I paw (x,y,z) 7 2 is described in overview. 

2 0 The raw data I raw (x,y,z) 7 2 is pre-processed 7 4 to produce 

a pre-processed angiographic volume I Q (x,y,z) 76 which 
conforms to a pre-selected data format and which has 
non-vascular contrast substantially removed therefrom. 
For example, the pre-processor 7 4 can convert the image 
25 format to one formed of isotropic voxels, pre-filter the 
data to remove noise, re-size the data, remove extraneous 
non-vascular contrast, and so forth. In the case of 
black blood angiographic (BBA) data, the pre-processor 74 
preferably removes non-vascular "black" contrast such as 

3 0 bone linings, air pockets, and muscle tissues that 

typically are imaged with low intensities similar to the 
black blood. The pre-processor 74 preferably also 
intensity- inverts BBA data so that the vascular regions 
correspond to high intensity areas of the pre-processed 
3 5 image I Q (x,y,z) 76. 

With continuing reference to FIGURE 5, the pre- 
processed volume I o (x,y,z) 76 is input to a three- 
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dimensional high curvature feature processor 78 which 
outputs a three-dimensional edge volume I edge (x,y,z) 80. 
The edge volume 80 emphasizes the edges of the vascular 
structures. Processing of angiographic data based on the 
5 edge volume 80 is advantageous because it retains the 
vessel lumen information and will accurately reflect 
narrowing of the blood vessels due to plaque buildup. 

With continuing reference to FIGURE 5 and with 
further reference to FIGURE 6A, a vessel centers processor 

10 82 receives the edge volume I eefee (x / y,z) 80 and processes it 
on a slice-by-slice basis to identify vessel centers 84 
substantially throughout the volume of interest which are 
representative of the imaged vascular system or systems. 
FIGURE 6A shows exemplary parallel slices 100 arranged 

15 perpendicular to the paper, with identified vessel centers 
marked or tagged by filled black circles 102. It will be 
appreciated that the shown vessel tags 102 are 
representative of a vessel portion having a lower portion 
104, a bifurcation point 106, and two vessel portions 108, 

20 110 extending upward from the bifurcation point 106. 
Although not immediately apparent from the two-dimensional 
schematic representation of FIGURE 6A, it will be 
appreciated that each parallel slice 100 extends two- 
dimensionally into and out of the paper and typically 

25 contains vessel centers across the two-dimensional 
surface, so that the combined slices 100 provide a set of 
vessel center tags 102 which are representative of the 
complete three-dimensional vascular system. 

With continuing reference to FIGURES 5 and 6A, 

30 the vessel center tags 84, 102 preferably include resolved 
vessel centers of overlapping vascular structures. In the 
exemplary schematic FIGURE 6A, the slice 100 o just above 
the bifurcation 106 includes two vessel centers 102,, 102 2 
which are likely to have overlapping vessel boundaries in 

3 5 the angiographic image. The vessel centers processor 82 
advantageously recognizes such overlapping conditions. 

The vessel center tags 84 form a valuable 
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database of global information about the vascular 
structures contained in the image. The tags 84 provide a 
valuable input for many contemplated angiographic post- 
acquisition processing applications. In one exemplary 
5 application (not shown) , the centers 84 are plotted to 
provide a graphical representation of the vascular system. 
Given a sufficiently high density of vessel centers 
corresponding to a close spacing of the parallel planes 
100, a reasonable representation of the vascular paths is 

10 provided thereby. This paths representation does not, 
however, include the vessel lumen information which is 
critical for identifying clinically significant plaque 
buildup and resultant vessel flow constriction which is 
typically of great interest to medical personnel. Neither 

15 does this simple approach provide for artery-vein 
separation or other sophisticated vascular analyses. 

A suitable processing method which incorporates 
the global vessel center tags 84 to effectuate rapid and 
accurate vessel segmentation including accurate vessel 

20 lumen reconstruction (even in the case of overlapping 
vessel images) and optional artery-vein separation is 
described with continuing reference to FIGURES 5 and with 
further reference to FIGURES 6B and 6C. A vessel 
segmentation processor 86 receives the tags 84. The 

25 processor 86 selects a starting location (x,y,z) 88 in the 
three-dimensional vascular system. The seed 88 can 
correspond to one of the tags 84. A direction processor 
90 receives the location 88 and analyzes the three- 
dimensional angiographic volume I Q (x,y,z) 76 in the 

30 neighborhood of the location 88 to determine the vessel 
direction 92 and the plane orthogonal to the vessel 
direction 94 at the location 88. The vessel segmentation 
processor 86 analyzes the orthogonal plane 94 to determine 
the vessel boundaries 120 (FIGURE 6B) in the plane 94, 

3 5 taking into account tagged vessel overlaps as necessary, 
and also calculates the next vessel center nearest to the 
point 88. The segmentation process is repeated until the 
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global set of vessel tags 84 are fully visited. With the 
vessel boundaries 120 determined, a three dimensional 
segmented vasculature 9 6 can be constructed, e.g. by- 
interpolation between the planes 124 and along the vessel 
5 centers 12 6 as shown in FIGURE 6C. 

By generating and employing the vessel center 
tags 84, the robustness of the vascular tracking system 7 0 
is greatly improved compared with trackers of the prior 
art. Prior art tracking systems operate in local space, 

10 working from a starting seed and tracking the direction 
locally. Such systems can fail at bifurcations, vessel 
gaps, vessel overlaps, and at other localized image 
structures which cause the vessel direction to become 
indeterminate. By incorporating the tags 84, a global 

15 database of information about the vasculature is included 
into the tracking. The tags 84 overcome the localized 
image structures by providing a global representation of 
the entire vascular system. For example, a vessel gap is 
identifiable by the existence of additional vessel center 

2 0 tags 84 positioned beyond the gap. Furthermore, the 
vessel tags identify overlapping vessels so that the 
tracker can account for the plurality of vessels at a 
point of overlapping or furcation. 

Having provided an overview of the tracking 

25 system embodiment 7 0 formed in accordance with the 
invention with reference to FIGURES 5 to 6C, the details 
of the exemplary segmentation system 7 0 are now described. 

With reference to FIGURE 7, a suitable 
embodiment of the pre-processor 7 4 is described. The raw 

30 angiographic volume I raw (x,y,z) 72 is input to an isotropic 
volume processor 140 which converts the gray scale volume 
image data 72, which may or may not be formed of cubic 
voxels, into an isotropic gray scale angiographic image 
volume 142 which is formed of cubic voxels of a selected 

35 resolution 144. Of course, if the raw data 72 as-acquired 
is isotropic, the conversion 140 is not applicable. 

With continuing reference to FIGURE 7, the 


-an- 
isotropic data is optionally down-sized using an 
edge-preserving processor 146 which in a suitable 
embodiment employs a shrinking and translating transform 
such as a discrete wavelet transform (DWT) 148 for edge- 
preserving re-sizing, to produce a re-sized image 150 
which is down-sized by a factor K 152. The DWT 
advantageously substantially preserves the high frequency 
edges of the image so that accurate blood vessel lumen 
information is retained. The wavelet transform is applied 
according to: 

<x> = 2~h(^) (1) 

where $ is the isotropic image volume 142 , a is a scale 
factor and b is a location/time parameter. The continuous 
wavelet transform is given by convolution as: 

-i. 

fja,b) = jf (t) 2~ 2 cs(!±) (2 ). 

— oo 

The mother function $ (herein the isotropic image volume 
142) should satisfy admissibility criterion to ensure that 
it is a localized zero-mean function. Equation (2) can be 
discretized by restraining a and b to a lattice of the 
form a=2" s to get a wavelet basis in which the dilations 
and translations (scaling and shifting) of the mother 
function $ is given by: 

s 

* {Sfl) (x) = 2 ? <3> (2" s x-i) < 3 > 

where s and 1 are integers that scale and dilate the 
mother function § to generate wavelets, such as those of 
the Daubechies wavelet family. The scale index s 
indicates the wavelet's width and the location index 1 
gives its position. Those skilled in the art will 
recognize that by applying equation (3) to the isotropic 
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image volume 142, the image is rescaled or dilated by K 
152 which is a power of 2, and translated by an integer. 
The image is transformed into four down-scaled images, 
three of which contain horizontal, diagonal, and vertical 
5 details, respectively, while the fourth image contains the 
desired downsized angiographic volume 150. 

In the case 154 of a black blood angiographic 
(BBA) volume, the BBA image is advantageously transformed 
by a pseudo-white-blood-angiographic (pseudo-WBA) 

10 processor 156 which: (i) replaces extraneous black non- 
vascular contrast with an alternative intensity which does 
not closely correspond with the black blood intensity; and 
(2) intensity-inverts the image so that the BBA volume has 
a "white-blood-type" contrast in which the vascular areas 

15 appear at high intensity. The processor 156 produces a 
pseudo-WBA image volume 158 which has the intensity 
polarity of WBA and is suitable for further processing in 
the same manner as a WBA image volume. 

With continuing reference to FIGURE 7, the image 

2 0 volume 150 (for WBA) or the image volume 158 (for BBA) is 

input to an intensity-stretching processor 160 which 
produces the pre-processed angiographic volume I o (x,y,z) 
76. This processor increases the intensity range by 
increasing the difference between the vessels intensity 
25 and the background intensity. In a suitable embodiment, 
the intensity-stretching processor 160 normalizes the 
voxel intensities according to: 

2 JV ~ 1 

^output ^ input * ~ J (4) 

max min 

where N is the number of bits per voxel, 1^ is the maximum 
intensity in a slice and I m1n is the minimum intensity in 

3 0 a slice. The intensity-stretching processor 160 improves 

the sensitivity to the vasculature by raising the vessels 
intensity above the background. It is most effectively 
used in conjunction with noise-filtering processors which 
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remove isolated, "hanging" noise voxels and otherwise 
suppress non-vascular intensities. 

With reference to FIGURES 8, a suitable 
embodiment of the pseudo-WBA processor 156 is described. 
5 The BBA volume I BBA (x,y,z) 150 BBA is processed to remove any 
non-vascular dark or black regions by constructing 17 0 a 
binary bone-air-muscle mask 172 and assigning 17 4 a gray 
intensity level to the non-vascular black areas identified 
by the mask 172. The gray intensity level assigned by the 

10 assignment processor 174 advantageously corresponds to the 
average value of the gray non-vascular areas which are not 
masked by the mask 172. Thus, the black non-vascular 
regions are replaced by a uniform average gray intensity. 

The output of the assignment processor 17 4 is a 

15 vessels-tissue gray scale BBA volume I vt (x,y,z) 17 6 which 
has clearly distinguishable vascular regions of very low 
intensity, i.e. black blood, on a higher intensity, gray 
background. The vessels-tissue gray scale BBA volume 
I vt (x,y,z) 176 advantageously replaces the non-vascular 

2 0 tissues which typically appear black in BBA images by a 
higher intensity corresponding to the average intensity of 
the gray areas of the BBA image. Those skilled in the art 
will recognize that the volume I vt (x,y,z) 176 overcomes the 
well-known difficulty in BBA image interpretation that 

2 5 certain non-vascular tissues, particularly air pockets, 

bone linings, and muscle tissues, appear black in BBA 
images and can interfere with, obscure, and confuse 
interpretation of the vascular information contained 
within the BBA image 150 BBA . 

3 0 With continuing reference to FIGURE 8, the 

vessels-tissue gray scale BBA volume I vt (x,y,z) 176 is 
intensity-inverted 178 to produce the pseudo-WBA image 
volume 158 in which blood vessels are represented by high 
intensity, i.e. appear white. With the intensity of a 
3 5 voxel of the gray scale BBA volume I vt (x,y,z) 17 6 
designated x vt and the corresponding voxel of the 
pseudo-WBA image volume 158 designated x yBA , the intensity 
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inversion is according to: 

x mA = (2 n -D -x vt (5) 

where n is the number of bits per voxel. For example, if 
n=8 then x WBA =2 55-x vt , and so a low intensity black pixel 
having x yt =0 converts to a high intensity x UBA =255. 
5 Although a linear voxel intensity inversion transformation 
has been described with reference to equation (5) , other 
intensity inversion transformations are also contemplated, 
such as non-linear transformations which advantageously 
account for non-linearities in the intensity data. 

10 With reference to FIGURE 9, a suitable 

embodiment of the bone-air-muscle mask processor 17 0 is 
described. A slice selector 190 selects a slice 192 from 
the BBA volume 150 BBA . Each pixel in the slice 192 is 
classified 194. The number of classes 196, designated 

15 herein as K, corresponds to the number of distinct types 
of pixels in the slice 192. As mentioned previously, 
typical BBA images contain two types of pixels: (1) black 
pixels that correspond to vascular structures or to 
certain non-vascular structures such as air pockets, bone 

2 0 linings, and muscle tissues; and (2) gray pixels that 
correspond to most organs and other tissue types and form 
much of the image background. Thus, the number of classes 
196 is selected as K=2 . 

With continuing reference to FIGURE 9, the class 

2 5 assignment processor 194, a suitable embodiment of which 
will be described later, assigns each pixel of the slice 
192 as either a black pixel or a gray pixel, corresponding 
to bone/air/vascular structures and tissue background 
respectively. The classification results are suitably 

30 expressed as an intermediate binary slice mask 198. The 
mask 198 is intermediate because it does not differentiate 
between vascular and non-vascular black pixels, whereas 
the desired mask should identify only the non-vascular 
black regions. 
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It will be appreciated that the vascular regions 
are tubular in structure, whereas non-vascular regions are 
typically large annular structures corresponding to bone 
linings or large irregular structures corresponding to air 
5 pockets. Thus, the mask processor 170 removes the 
vascular regions from the intermediate mask 198 using a 
mathematical reduction/ enlargement seguence 200. In a 
suitable embodiment, an erosion process applies a moving 
window 202 to erode the edges of the black regions of the 

10 mask by a pre-selected amount. The narrow tubular 
vascular structures are thus completely removed whereas 
the non-vascular structures, which are typically larger, 
have only the edges reduced. A subsequent dilation 
process enlarges the non-vascular areas back to their 

15 original uneroded size. The resulting slice mask 204 
contains only the non-vascular structures with the blood 
vessels removed, and is suitable for use in identifying 
the non-vascular black regions of the image. 

With continuing reference to FIGURE 9, the slice 

2 0 mask generation process is repeated on a slice-by-slice 
basis 206, 208 to produce a mask slice 204 corresponding 
to each slice of the BBA volume 150^. The mask slices 204 
are combined 210 to generate a mask volume 172. 

With reference to FIGURE 10, a suitable 

25 embodiment of the class assignment processor 194 is 
described. The described processor 194 embodiment is 
based on constructing a parameterized statistical model 
including a Gaussian distribution. The mean and standard 
deviation of the Gaussian distribution are the optimized 

30 parameters. The class assignment processor 194 is 
described with reference to N pixels in the image slice 
192, and with respect to the K classes 196. Although K=2 
for typical two-class BBA data, inclusion of additional 
classes is also contemplated for specific imaging 

35 situations, and so the class assignment processor 194 is 
described for an arbitrary number of classes K 196. The 
pixels are designated herein by i where l<i<N, and the 
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classes are designated by k where l<k<K. The pixel 
intensities are expressed herein as x 1 . The pixel 
classifications are expressed as y 1 - k , which indicates the 
statistical probability that the pixel i is a member of 
the class k. A Gaussian probability model 220 is used to 
calculate unweighted pixel probabilities 222 according to: 


Pi, k =Plx 1 \y itk > *) = . e 2o > (6) 

/(2n) D o^ 


where: P ifk =P (*,- [ y,- , $) is the unweighted Bayesian 
probability that a pixel of intensity x i is in the class k; 
M k and a k are the mean and standard deviation of the 

10 Gaussian probability distribution model; $ are the model 
parameters, e.g. *={/V a 0 ' and D is a dimensional 
parameter, here D=l. 

With continuing reference to FIGURE 10, the 
described embodiment of the class assignment processor 194 

15 optimizes the classification probabilities by maximizing 
the a posteriori probability by iteratively optimizing the 
Gaussian model parameters §={^ k ,a k ) and a weighting w fc . The 
pixel classification values y. k and the distribution 
parameters $={^ k ,a k } are assigned initial values 224, and 

2 0 an initial probability weighting w k is computed. Suitable 
initial values are given by: 

w°= 1 F o_ (k+1 ) 2 N 0 _ In n _ n lns 


where E, is the covariance and the superscript "o" indicates 
the zeroth iteration. The initial mean and standard 
deviation 224 are copied 226 into the past-iteration or 
25 old mean and standard deviation memory 228. Using the 
Gaussian probability model 220, the unweighted pixel 
probability p ik 222 is calculated using the old mean and 
standard deviation 228, and an updating of the weighted 


I* 
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n+1 


classification values y,- k is performed 230 according to: 

n n 
n+l_ Pi,k W k 

y ^ fc ~^ ^ n <8), 


where the superscript n indicates the current iteration, 
the superscript n+1 indicates the next iteration, and w k n 
is a weighting given by: 

En-l 
Yi,k 


n 


W k = ( 9 ) 


5 The updated classification probabilities y. . n+1 232, which 

1 t K 

are weighted probabilities, are thus obtained. An 
iterative updating of the statistical model parameters is 
also performed 234 according to: 


vT - 


i 

H Yi,k 


n+1 


^ 


i 

12 y/!* 


(10) 


to generate updated statistical model parameters 23 6. 

10 With continuing reference to FIGURE 10, the 

iterative process is terminated using a suitable 
convergence criteria. For example, convergence can be 
indicated by computing 2 38 a difference between the old 
mean /it k n and the new mean /i k ^ 1 for each class k and checking 

15 240 for convergence by comparing the differences with pre- 
set criteria. The iterative process can also optionally 
be terminated after a pre-set number of iterations. 

Once the iterative process is complete, an 
assignment optimizer 242 assigns the optimized 

20 classification to each pixel of the intermediate slice 
mask 198. For each pixel, the final weighted 
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probabilities y, v of each class for a pixel are obtained 
and the class k that has the highest corresponding 
weighted probability y, k for that pixel is selected as the 
classification of that pixel. 
5 A suitable embodiment of an exemplary 

pre-processor 7 4 has been described with reference to 
FIGURES 7-10. The exemplary pre-processor 7 4 generates an 
isotropic, edge-preserving, optionally re-sized and 
intensity- inverted image I Q (x,y,z) 76 which has WBA 

10 intensity polarity. Of course, additional or other pre- 
processing steps are also contemplated, such as additional 
noise filtering, smoothing, etc. 

With reference returning to FIGURE 5, the pre- 
processed angiographic volume I Q (x,y,z) 76 is input to the 

15 three-dimensional high curvature feature processor 7 8 
which outputs an edge volume I edge (x,y, z) 80. 

With reference to FIGURE 11, a suitable 
embodiment of the edge processor 7 8 is described. The 
edge processor 78 approximates differentiation of the 

20 angiographic volume I o (x,y,z) 76 by convolving the image 76 
with a second order Gaussian derivative. A voxel 2 62 is 
chosen 2 60. A Gaussian design processor 2 64 calculates a 
Gaussian and its derivatives 2 66 based on initial 
coefficients 268. The calculation of the Gaussian 

25 derivatives preferably takes advantage of known closed- 
form expressions for the Gaussian derivatives. Although 
the convolving is ordinarily a calculation of order 0(k 3 ) 
where k is the number of number of voxels forming the 
Gaussian kernel 2 66, by convolving in separable mode the 

3 0 calculation is reduced to order 0(3k) . Thus, the 
convolving reduces to: a convolution in the x-direction 
270 which produces an x- intermediate result 272; a 
convolution in the y-direction 27 4 which produces a y- 
intermediate result 27 6; and a convolution in the z- 

35 direction 278 which produces the final result 280. The 
gradient processor 282 combines the one-dimensional 
derivatives estimated by the convolving to form an edge 
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volume voxel according to: 


J edge (X ' y ' Z > 









( ay, 


{ dzj 


(11) 


where I is the pre-processed angiographic volume I o (x,y,z) 
7 6 and the partial derivatives are evaluated at the 
location (x,y,z) of the selected voxel 262. After 
5 repeating 284 for every voxel, the edge volume I edge ( x /y/ z ) 
80 is generated. 

It will be appreciated that the Gaussian is a 
function with a variance a. If a is large, the Gaussian 
convolution can take prohibitively long, especially for 

10 higher order derivatives. For calculation of the edge 
volume 80, second order derivatives are used. In a 
preferred embodiment which greatly improves calculation 
speed, the convolving is performed using a z-transform 
decomposition of the Gaussian into sinusoidal and 

15 cosinusoidal components. 

With reference to FIGURE 12, a suitable 
embodiment for performing the z-transform decomposition is 
described. For a scale a 290 and a one-dimensional 
coordinate x 292, a cosine processor 294 calculates a 

20 cosine term 296, and similarly a sine processor 298 
calculates a sine term 300. The cosine and sine terms 
are multiplied 302, 304 by the cosine and sine z-transform 
coefficients 306, 308, respectively, to produce the cosine 
and sine z-transform components 310, 312, respectively. 

2 5 These components 310, 312 are additively combined 314 to 

produce the z-transform 316. 

With reference to FIGURE 13, an exemplary time- 
of-flight white blood angiographic image 330 of the 
carotid area is shown. The image 330 includes five 

3 0 isolated vessels 332 as well as one pair of overlapping 

vessels 334. 

With reference to FIGURE 14, the edge volume 340 
corresponding to the image 330 is shown. The edge volume 
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is formed according to the method described with reference 
to FIGURES 11 and 12. Each of the isolated vessels 332 is 
shown to have a well defined closed contour 342 
corresponding thereto. Similarly, the overlapping vessels 
5 pair 334 has a corresponding well defined closed contour 
344. Besides the contours 342, 344, the remainder of the 
edge volume 340 is substantially free of features, except 
for a few low intensity enclosed regions 346 associated 
with non-vascular features. These non-vascular edges 346 

10 are lower intensity versus the vessel contours 342, 344. 

With reference back to FIGURE 5, the vessel 
centers processor 82 receives the edge volume 80 and 
identifies a plurality of vessel center tags 84 
representative of the vascular system. The vessel center 

15 tags 84 identify the vessel centers in three-dimensional 
space, e.g. by coordinates (x,y,z), and also 
advantageously identify vessel centers corresponding to 
overlapping vessel images. 

With reference to FIGURE 15, a suitable 

2 0 embodiment of the vessel centers processor 82 is 

described. The approach of this embodiment is to divide 
37 0 the three-dimensional edge volume 80 into two- 
dimensional edge volume slices 372. Each slice 376 is 
selected 37 4 for processing in turn. The edge volume 
25 slice 376 is converted 378 into a binary slice 380 by 
setting the non-zero values to one 382. The binarization 
378 can be effectuated, for example, by using a very low 
thresholding value. As seen in the exemplary edge volume 
slice 340 of FIGURE 14, the edge volume image elements are 

3 0 substantially zero (i.e., black) except at the edge 

contours 342, 344, and so a low thresholding value will 
appropriately set the vessel contours 342, 344 to binary 
one imposed on a binary zero background, thus forming 
closed binary contour edges of the vessel structures. The 
35 vessel structures are flood-filled 384 to form a slice 
having flood-filled vessel structures 386. There are a 
number of well-known methods for flood-filling closed 
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contours, and one or more appropriate algorithms are 
typically available as standard functions in signal 
processing software development environments. The flood- 
filled regions are counted 388 to produce a count and 
5 location 390 of the vessel structures in the slice 386. 
In a suitable embodiment, a run length encoder method 
performs the counting. Run length encoders are well-known 
algorithms which are commonly employed in JPEG image 
compression, and functions which implement a run length 

10 encoder or a variation thereof are often available in 
signal processing software toolboxes or libraries. 

With continuing reference to FIGURE 15 and with 
further reference to FIGURE 16, a vessel centers tagger 
400 receives the slice 386 with the flood filled vessel 

15 structures and the count and location of the vessel 
structures 390. For each flood-filled vessel structure of 
the slice 386, the vessel centers tagger 400 identifies 
one or more vessel centers 402 corresponding thereto. The 
vessel centers tagger 400 is suitably embodied by a 

20 recursive erosion process to be described later. 
Advantageously, the vessel centers processor 82 also 
particularly tags overlapping vessel structures by 
identifying 404 vessel structures having more than one 
vessel center 402, since a plurality of vessel centers 402 

25 corresponding to a single vessel structure indicates the 
singular vessel structure is formed of overlapping vessel 
images. Appropriate vessel overlap tags 406 are generated 
by the vessel overlap identifier 404. By cycling 408 
through the slices 372 a set of vessel centers and overlap 

3 0 tags 84 is generated that is representative of the imaged 
vascular system. 

With reference to FIGURES 17A to 17D, a suitable 
method for implementing the vessel centers tagger 400 is 
described. A flood filled vessel structure 420 is 

35 successively eroded as shown in FIGURES 17B through 17D 
using a moving window 422 . The recursive erosion process 
iteratively removes outer portions 424 of the vessel 


- 30 - 

structure 420, leaving successively smaller vessel 
structure portions 426, 428. The final remaining vessel 
structure portion 430 after the recursive eroding tagged 
as the vessel center, e.g. its (x,y,z) spatial coordinates 
5 are recorded in the table of vessel centers 84. It will 
be recognized that the flood filled vessel structure 420 
of FIGURE 17A erodes to a single vessel center 430. 

With reference to FIGURES 18A to 18 D, the method 
of FIGURES 17A to 17D is applied to a flood filled vessel 

10 structure 440 which images overlapping vessels. During 
the initial stages of the recursive eroding, the moving 
window 422 removes outer layers 444 leaving a single 
reduced vessel structure portion 446. However, as the 
iterative eroding continues, two distinct vessel structure 

15 portions 448, 450 become separately distinguishable. 
Ultimately, the vessel structure portions 448, 450 erode 
down to identifiable vessel centers 452, 454. Since the 
vessel centers 452, 454 correspond to a single flood 
filled vessel structure 440, the vessel centers processor 

20 82 identifies 404 the flood-filled image 440 formed of 
overlapping vessel images corresponding to the vessel 
centers 452, 454. 

The moving window 422 shown in exemplary FIGURES 
17B, 17C, 18B, and 18C is a filled circular element. With 

25 reference to FIGURE 19, an exemplary square window 470 is 
shown which is formed from nine pixel elements arranged 
about a center (r,c). In order to improve the speed of 
the recursive erosion, the erosion is advantageously 
separated into two passes. A first pass runs across the 

3 0 slice from a first corner, e.g. the upper left corner, to 
a second opposite corner, e.g. the lower right corner, and 
uses only the upper left portion 472 of the moving window 
470. The first pass is followed by a second pass that 
runs across the slice from the second corner, e.g. the 

3 5 lower right corner, to the first corner, e.g. the upper 
left corner, and uses only the lower right portion 47 4 of 
the moving window 470. 
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With reference to FIGURE 20, the first pass of 
the decomposed recursive erosion embodiment of the vessel 
centers tagger 400 is described. For each pixel (r,c) of 
the slice 386 having flood-filled vessel structures, the 
5 top four values A, B, C, D around (r,c) are obtained 500, 
namely: A=(r-1, c-1) ; B=(r-1, c) ; C=(r-1, c+i) ; and D=(r, 
c-1) . It will be appreciated that these top four values 
502 correspond to the upper left portion 472 of the moving 
window 470 of FIGURE 19. A minimum processor 504 selects 
10 the minimum value 50 6 of these four values. And 
add/replacement processor 508 adds one to the minimum 
value 506 and replaces the pixel (r,c) with the new value 
510. Thus, the transform for the first pass is: 

(r,c) new = rain {A, B,C,D} + 1 (i 2 ) 

where A, B, C, D are the pixels defined above. The 
15 transform of equation (12) is repeated 512, 514 for each 
pixel to produce an intermediate LRTB image 516. 

With reference to FIGURE 21, the second pass of 
the decomposed recursive erosion embodiment of the vessel 
centers tagger 400 is described. For each pixel (r,c) of 
20 the intermediate LRTB image 516, the bottom five values 
around (r,c) are obtained 530, namely: E=(r, c) ; F=(r+1, 
c+1) ; G=(r, c+1) ; H=(r+1, c) ; and I=(r+1, c-1). It will 
be appreciated that these bottom five values 532 
correspond to the lower right portion 47 4 of the moving 
25 window 470 of FIGURE 19. An add/replacement processor 534 
adds one to each of the five values E, F, G, H, I to 
create five new values 53 6. A minimum processor 538 
selects the minimum value as the new value 54 0 at the 
pixel (r,c). Thus, the transform for the second pass is: 

(r,c) new = m±n {E+l , F+l ,G+1 ,H+1 , 1+1} ( 13) 

30 where E, F, G, H, I are the pixels defined above. The 
transform of equation (13) is repeated 542, 544 for each 
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pixel to produce the eroded image 546 resulting from an 
erosion pass. The recursive erosion process of FIGURES 20 
and 21 is repeated until the final vessel centers are 
obtained . 

5 With reference to FIGURES 2 2A through 2 2D, 

exemplary results for recursive erosion according to the 
process of FIGURES 20 and 21 are presented. In FIGURE 
21A, a synthetic structure 57 0 is formed of two fully 
distinct flood-filled circular contours 572, 574. As 

10 shown in FIGURE 22B, the recursive erosion process 
performed according to FIGURES 2 0 and 21 produces two 
eroded structures 576, 578 that have clearly defined high- 
intensity centers corresponding to the centers of the two 
flood-filled circular contours 572, 574. 

15 With reference to FIGURES 22C through 22D, a 

synthetic structure 580 is formed of two fully distinct 
flood-filled circular contours 582, 584. As shown in 
FIGURE 22D, the recursive erosion process performed 
according to FIGURES 2 0 and 21 produces two eroded 

2 0 structures 586, 588 that have clearly defined high- 
intensity centers corresponding to the centers of the two 
flood-filled circular contours 582, 584. A high intensity 
connecting line 590 connects the high intensity centers 
58 6, 588 but does not detract from the detectability of 

2 5 the centers 586, 588 because they are positioned at the 

endpoints of the line 590. 

With reference to FIGURE 23, the exemplary 
vessel segmentation processor described in brief overview 
with reference to FIGURE 5 is described in greater detail. 

3 0 The segmentation process of FIGURE 23 uses a path tracking 

process 600. A starting point 610 lying within the 
vascular system is selected. The starting point is 
typically a root of the venous or arterial system of 
interest, and can optionally be selected manually or 
35 identified using an automated system, e.g. by selecting 
the vessel structure having the largest area from the 
count and location of vessel structures 390 (FIGURE 15) . 
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The starting point is optionally selected from the table 
of vessel center tags 84. The tracking point at (x,y,z) 
88 is selected 612, initially corresponding to the 
starting point 610. 
5 Tracking from the point (x,y,z) 88 includes (1) 

finding the vessel direction; (2) finding the vessel 
center; (3) finding the extent of the vessel about the 
vessel center in the plane orthogonal to the vessel 
direction; (4) locating the next vessel center; and (5) 

10 repeating the process. The direction processor 90 finds 
the vessel direction 92 and the orthogonal plane 94. A 
vessel centers search processor 612 searches the vessel 
centers tags 84 to locate the nearest vessel center 614 in 
space to the point (x,y,z) 88. (Of course, if the point 

15 88 is selected from the vessel centers tags 84, the 
searching 612 can be omitted) . The segmentation processor 
8 6 analyzes the orthogonal plane 94 to determine the 
extent of the vessel about the vessel center 614. The 
boundaries of the vessel along with the vessel center 612 

20 form a vessel skeleton portion 616. The process is 
repeated 618, by selecting 612 a new point 614. The 
selecting 612 advantageously moves a pre-selected distance 
along the vessel direction 92, e.g. a distance equivalent 
to one voxel for the best tracking resolution. Once the 

25 vessel termination is reached, a vascular path 96' is 
obtained. 

The tracking system described with reference to 
FIGURE 23 applies to a vascular path 96'. At bifurcation 
points, the overlap tags 406 (FIGURE 16) are 

3 0 advantageously used to identify new starting points 610 at 
which the path tracking 600 of FIGURE 23 is repeated. 

With reference to FIGURE 24, a robust general 
tracking system 630 that uses the tags 84 to ensure 
tracking of all the vascular branches is described. The 

3 5 vessel starting point 610 is selected 632 from the vessel 
center tags 84. The tracking system 600 of FIGURE 2 3 is 
applied to generate the tracked path 96' corresponding to 
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the starting point 610. During the tracking 600, each 
tracked vessel center 614 is marked as visited. After the 
tracking 610, a check 634 is performed to determine 
whether all the vessel tags 84 have been marked as 
5 visited. If there are unvisited tags, one of the 
unvisited tags is selected 632 as the new starting point 
610, and the tracking 600 is repeated. Once all the 
vessel centers 84 have been visited, the tracked paths 96' 
are combined to form the full vascular skeleton 96. 

10 Because the vessel center tags 84 are representative of 
the imaged vascular system, the system is robust. 
Bifurcation points, tortuous or occluded vessels, vessel 
overlaps, intertwined vessels, partial volume averaging 
and other imaging artifacts, and vessel gaps which act 

15 alone or in various combinations to produce localized 
disjoints or multiplicities of the vascular path do not 
prevent the general tracking system 630 from tracking all 
vessel branches, because visitation of every vessel branch 
is ensured by visiting every one of the vessel tags 84. 

2 0 With reference to FIGURE 25, a suitable 

embodiment of the direction processor 90 is described. 
The concepts of differential geometry, curves and surfaces 
are preferably used for estimating the direction 92 and 
the orthogonal plane 94. In 2-D images, lines or curves 
25 have two directions, one direction is along the line or 
curve and second direction that is orthogonal to the line 
or curve. These directions are the vectors corresponding 
to the eigenvalues of the Hessian matrix at a given point. 
For a 2-D image, the two eigenvalues are the two 

3 0 curvatures at that point, i.e. the maximum curvature (k,) 

and the minimum curvature (k 2 ) . Using the differential 
geometry, one can represent the amplitude of the curvature 
and its corresponding direction by the eigenvalues and 
eigenvectors of the Weingarten matrix W = F 1 " 1 F 2 , where F 1 
3 5 is a function of the partial derivatives of the image in 
x and y directions, and F 2 is the function of the second 
partial derivatives in x and y directions. 
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The three dimensional pre-processed angiographic 
volume I o (x,y,z) 76 is modeled as a four-dimensional 
surface given by: 


S = 


x 

y 

z 

I(x,y, z) 


(14) 


10 


where the first three elements are the spatial coordinates 
(x,y,z) and the fourth element is the image intensity 
I(x,y,z) at the voxel location (x,y,z). In this case, 
there are three principal curvatures. Two curvatures 
correspond to the two orthogonal directions in the cross- 
sectional plane of the tubular blood vessel, while the 
third principal curvature corresponds to the vessel 
direction or vessel orientation. The three directions can 
be computed using the Weingarten matrix which is a 
combination of the first and second form of the 


hypersurf ace, i.e. W = F. 


2' 


where F. is the fundamental 


15 form of hypersurf ace and a function of I. I and I , the 

x y z 

three partial derivatives of the image volume, and F 2 is 
the second fundamental form of hypersurface and is a 
combination of the second partial derivatives I xx , 


I xy' "'"xz' 


I yy / I yz / and I zz . More explicitly, 


W = F X X F 2 


where 


F, = 


1+Iy, II II 

x x y x z 


ii 1 + i i 

x y ¥ y z 


II II 1+1' 

x z y z 2 


F„ = 


III 

xx xy xz 

III 

xy yy yz 

III 

xz yz zz 


(15) 


2 0 A Weingarten (W) matrix 652 is constructed 650 according 
to equation (15) . The eigenvectors and eigenvalues 656 of 
this matrix are obtained by diagonalizing 654 the 
Weingarten matrix, such eigenvalue decomposition being 
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well known to those of ordinary skill in the art. The 
vessel direction 92 and the orthogonal plane 9 4 are 
identified 658 from the eigenvalues and eigenvectors 656 
as follows. The eigenvector direction corresponding to 
5 the minimum curvature is the vessel direction 92 . The 
plane passing though the point (x,y,z) whose normal is the 
vessel direction 92 is the orthogonal plane 94. 

The orthogonal plane 94 is in general an obligue 
plane relative to the principle axes of the isotropic MRA 

10 volume 76. That is, the plane 94 is not in general 
parallel to one of the axial, sagittal, and coronal planes 
along which MRA data is typically measured. Trilinear 
interpolation of the MRA volume 76 is preferably performed 
to obtain the orthogonal plane 94 with isotropic pixels. 

15 With reference returning to FIGURE 23, the 

vessel segmentation processor 86 receives the orthogonal 
plane 94 and estimates the vessel boundaries in the 
orthogonal plane 94 which form part of the vessel skeleton 
616. 

20 With reference to FIGURE 26A and 26B, a suitable 

embodiment of a method for estimating the vessel 
boundaries is schematically described. FIGURES 2 6A and 
26B schematically show estimation of the vessel boundaries 
for the flood-filled overlapping vessel structure 440 of 

25 FIGURE 18A. It is known from the table of vessel centers 
and overlap tags 84 that there are two vessel centers 452, 
4 54 associated with the overlapping vessel structure 440 
(see FIGURE 17D) . Initial geometric contours 680, 682 are 
arranged around the vessel centers as shown in FIGURE 26A. 

30 The contours 680, 682 are level set contours defined by 
the partial differential equation: 

|£ = (sk + V p ) |Vct>| - V ext Vd> (16) 


where <p is the level set contour 680, 682, ck, V p , and V ext 
are speed functions known as the curvature, regional, and 
gradient speed functions, respectively, that control the 
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propagation or flow of the contour. The level set 
contours 680, 682 are iteratively expanded from their 
initial position based on the speed functions according to 
equation (16) , and as shown in FIGURE 2 6B, to produces 
5 expanded level set contours 684, 686 which iteratively 
fill and define the vessel boundaries. The region-based 
level set approach has been previously used in brain 
tissue segmentation applications. Its application herein 
to vascular boundary estimation provides improved accuracy 

10 and speed. Methods previously used for vessel boundary 
estimation, particularly parametric contour methods, are 
inaccurate at sharp corners and have a tendency to bleed 
through gaps in the vessel structure. 

With reference to FIGURE 27, the level set 

15 approach as applied to vascular boundary reconstruction is 
described. An initial contour 702 is constructed 7 00 
about the vessel center selected from the vessel center 
tags 84 in the orthogonal plane 94. The vessel contour is 
propagated, i.e. expanded 704. The propagating 704 is 

2 0 iterative, occurs within the plane 94, and is constrained 

by appropriate constraints 706 such as the distance from 
the vessel center and contact between the propagating 7 04 
contour and neighboring contours. Another constraint is 
the outer boundary 708 of the vessel structure in the 
25 image. These constraints serve as stopping criteria which 
prevent the contour from expanding too far, e.g. beyond 
the vessel structure or into another vessel. The final 
propagated contour defines the vessel boundary 710. The 
level set approach is repeated 712 for each vessel center 

3 0 of the vascular structure. For the exemplary case of the 

vessel structure 440 (FIGURE 18A) the level set approach 
is repeated 712 for the two vessel centers 452, 454 as 
shown in FIGURES 2 6A and 2 6B. 

With reference to FIGURE 28, a suitable 
35 embodiment of the contour propagator 704 is described. 
The propagator 704 of FIGURE 28 is described with 
reference to a pair of overlapping vessels such as are 


A 
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shown in FIGURES 2 6A and 2 6B. Those skilled in the art 
can easily make the necessary changes to accommodate a 
vessel structure having only a single vessel center 
associated therewith. Thus, there are two initial 
5 contours 702, for which initial level set distance map 
fields 732 are computed 730 based on: the arrangement of 
the initial contours 7 02 about the vessel centers; the 
orientation of the orthogonal image plane 94; and the 
selected width of a narrow band 734 within which the level 
10 set field 732 is applied. 

With continuing reference to FIGURE 28, a new 
field 738 is calculated 736 by solving equation (16) 
iteratively to yield: 

<C = < y " At {V (x,tf +V grad (x,y) -^ cur (x,y) } (17) 


where an integrated speed input 740 includes the three 
15 velocity terms V reg (x,y), V grad (x,y), and V cur (x,y). 
Typically, At=l, i.e. each iteration corresponds to a 
discrete time interval. V peg (x,y) is a regional speed or 
velocity given by: 

V e? (x,y) = max{V p , 0}V + + min{V p , 0 } V" 
where 

co„ (18), 

V (x,y) = 

p '* Y [l-2u(x,y)] 

v* = [v>v;i 1/2 , v- = [v; + v;] 1/2 


V ,(x,y) is a gradient speed or velocity given by: 
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V w (x,y) = V . (x,y) + V . (x, y) 

gr<3d 1 ' J ' grad, x v ' ' grad, y 1 

(19), 

max{p n (x,y) ,0} £> " x ( x, y) +min{ g n ( x, y) , 0 } D +X (x,y) , 

V . {x,y) = 

grad, y ' J ' 

max {g n (x,y) , 0} D " y (x, y) +min{p n ( x,y) , 0} D +y (x,y) 
and V grad (x,y) is a curvature speed or velocity given by: 
v cur^ x 'y^ =eK"(x,y) [ (D 0x (x,y) ) 2 + (D 0y (x,y) ) 2 ] 1/2 (20). 

In the above equations, y is a damping coefficient, u(x,y) 
is a fuzzy pixel membership function whose value lies in 
the range [0,1], u R is a regional weighting, p n and q n are 
5 the x- and y-components of the gradient strength and are 
calculated from the edge volume 80, K n (x,y) is the 
curvature at location (x,y) , and V v + , V + , V", V " are the 

a y a y 

forward and backward level set gradients in the x and y 
directions which can be expressed in terms of finite 

10 difference operators of the form D +/ " <x,y> (x , y ) for forward 
and backward differences, and as D 0<x ' y> (x, y) for averaged 
differences. The regional velocity V reg relates to the 
fuzzy membership function u(x,y) which for example could 
be computed using the fuzzy C mean algorithm known to the 

15 art. The gradient velocity V grad relates to the gray scale 
intensity gradients, and the curvature velocity V cur relates 
to the curvature in level set space. 

With continuing reference to FIGURE 28, once the 
new level set field calculated according to equation (17) , 

2 0 the contour is updated and the level set field is 
reinitialized 742 using isocontouring 746 in which the new 
contour is defined essentially corresponding to the zero 
contour of the level set field 738. The updating can 
employ the fast marching method or other techniques known 

25 to the art. The calculating 742 of the updated contour 
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and field 744 is also constrained by selected constraints 
706 as discussed previously with reference to FIGURE 27. 
Of course, since FIGURE 28 refers to an exemplary 
propagation of a pair of overlapping vessels such as are 
5 shown in FIGURES 26A and 26B, there are two contours 684, 
68 6 (FIGURE 2 6B) being calculated, one each about the 
vessel centers 452, 454 (FIGURE 18D) . Thus, an 

appropriate constraint in this exemplary case is that the 
two contours 684, 686 not overlap. Similarly, if there 
10 were three or more contours corresponding to an 
overlapping vessel structure, an appropriate constraint 
would be the non-overlapping of the three-or more 
contours . 

With continuing reference to FIGURE 28, the 
15 iterative level set propagating is repeated 748, 750, 752 
a plurality of times until the contours converge to the 
two zero-level contours 710 corresponding to the two 
overlapping vessels. 

Although a contour propagation employing a 
2 0 geometric contour in region-based level set space 
including a regional fuzzy membership velocity component 
is described herein, other methods for finding the vessel 
boundaries are also contemplated, such as methods using 
parametric contours, geometric contours employing other 

2 5 types of distance maps, and the like. 

With reference returning to FIGURE 24, the 
exemplary robust vessel tracking system 630 generates a 
full vascular skeleton segmentation 96. This extracted 
information can be conveyed to medical personnel in a 

3 0 variety of ways, including three-dimensional graphical 

displays, two-dimensional projections, selected vascular 
sub-system displays (e.g., limited to the arteriogram), 
and the like. 

With continuing reference to FIGURE 2 4 and with 
35 further reference to FIGURE 29, an exemplary embodiment of 
a method for generating a suitable three-dimensional 
display is described. The vascular skeleton 9 6 is 
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supplied as an input. The tracking system 630 generates 
a skeleton 96 formed of a sequence of vessel contours or 
boundaries 710 (FIGURES 27 and 28) lying in two 
dimensional slices corresponding to the orthogonal planes 
5 94 (FIGURE 23) . Thus, the first step is to operatively 
interpolate across the boundaries 710. In the exemplary 
embodiment of FIGURE 29, shaded wire mesh method of a type 
known to the art is employed. A triangulation method 77 0 
generates a three-dimensional wire mesh of triangular 

10 elements 772. The normals of the vertex triangles 776 
which determine the solid vascular shape are computed 774. 
A display processor 778 converts the wire mesh into a 
shaded three-dimensional display 780. Optionally, 
selected graphical animation is produced 778 as well, such 

15 as graphical highlighting of a selected vascular system, 
e.g. the arteriogram. Of course, the three-dimensional 
display just described is exemplary only, and those 
skilled in the art will recognize that many other methods 
of conveying the segmented vascular tree information can 

20 also be used. 

The invention has been described with reference 
to the preferred embodiments. Obviously, modifications 
and alterations will occur to others upon reading and 
understanding the preceding detailed description. It is 

2 5 intended that the invention be construed as including all 
such modifications and alterations insofar as they come 
within the scope of the appended claims or the equivalents 
thereof . 


